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Abstract 



We numerically test the mode-coupling model (Lubow 1991a) of tidal 
instability in SU UMa systems. So far, all numercial models confirming it 
have been based on SPH codes and isothermal equation of state. In our 
paper we present Eulerians models, using both isothermal approximation 
and the full energy equation. We also investigate influence of different 
ways of mass transfer. While isothermal models behave similarly to SPH 
simulations, the behaviur of models with full energy equation is quite 
different, and the mode-coupling model is not confirmed in this case. 

1 Introduction 

SU UMa stars form a subclass of dwarf novae, which, in turn, are a subclass of 
cataclysmic variables (CVs). Like all CVs, the SU UMa stars are semidetached 
binary systems consisting of a white dwarf (the primary) and a low-mass, main- 
sequence star (the secondary). The secondary fills its Roche lobe, and a stream 
of gas flows from its surface toward the primary through the inner Lagrangian 
point. Because of excessive angular momentum, the stream is deflected from 
its original direction, and an accretion disk is formed around the primary. The 
disk may be subject to a thermal instability (Smak 1999), resulting in episodes 
of enhanced accretion rate. To a distant obsever, such an episode is visible as a 
temporary brightening of the star, commonly referred to as an outburst. 

As opposed to ordinary dwarf novae, the SU UMa stars exhibit a clearly 
bimodal distribution of outbursts. Normal outbursts have an amplitude of ~ 3 
mag, and last from one to four days. Superoutbursts are by ^ 1 mag stronger, 
and last for up to several weeks. The recurrence time of normal outbursts (days 
to weeks) is not constant, and it varies substantially form one system to another. 
The superoutbursts repeat more regularly, and their recurrence time is much 
longer (months to years). In the extreme case of WZ Sge stars, superoutbursts 
are nearly exclusively observed, with a recurrence time of up to several tens of 
years. In superoutbursts, the light curve of a SU UMa system is modulated with 
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a period a few per cent longer than the orbital period. Those modulations are 
referred to as super humps. The super hump signal is known to originate from 
extended source(s) in the outer disk (Warner 1995). 

Superoutbursts are thought to be driven by a combination of thermal and 
tidal instability. During normal outbursts, the disk grows in size as it diffuses 
under the influence of increased viscosity. In systems with mass ratios g < 0.25, 
it eventually reaches up to the location of the 3:1 resonance, at which the orbital 
frequency of the disk gas is three times larger than the orbital frequency of 
the binary (we define q as the ratio M2/M1, where Mi and M2 are primary's 
and secondary's mass, respectively). Subsequently, the tidal instability sets 
in, the disk becomes eccentric, and, seen in the inertial frame, it performs a 
slow, prograde precession. The tidal influence of the secondary on (and the 
viscous dissipation in) the outer disk, is largest when the bulk of the disk passes 
the secondary. The superhump period is then the beat period between the 
precession period and the orbital period of the binary (Osaki 1996). 

Disk precession and superhump phenomenon have been subject to rather 
intense theoretical investigations, largely based on numerical simulations. Ac- 
cording to Lubow (1991a, b), the eccentricity builds up due to nonlinear inter- 
action of waves, in which the m=3 component of the tidal field is a key factor. 
Heemskerk (1994) performed simulations using only that component, and he 
found that the disk became eccentric, but it precessed retrogradely. Moreover, 
with the full tidal potential, the accretion disk was kept away from the location 
of the resonance, and no significant eccentricity was produced. Heemskerk's 
results are the only ones obtained with an Eulerian (fluid) code. All remaining 
models presented in the literature have been based on Lagrangian (particle) 
codes. A detailed review of those calculations can be found in Murray (1998). 
While a qualitative agreement with observations of superhump systems was 
reached in several aspects, the models did not entirely agree with the analytical 
theory of the tidal instability. In particular, the measured eccentricity growth 
rates were much smaller than the predicted ones. The models themselves were 
often based on an extremely simple physical scenario (fully isothermal disk; gas 
from the secondary uniformly "raining" onto a circular orbit within the disk). 
In some of them the tidal instability was initiated by an arbitrary increase of 
viscosity in the disk by a factor of 10. Finally, even in the models which properly 
followed the stream of gas between the inner Lagrangian point up to its collision 
with the edge of the disk, the resolution in the collision region was too poor to 
resolve strong shock waves responsible for the hot spot phenomenon. 

In the present paper, we obtain Eulerian models of disks in SU Uma systems 
in order to isolate the influence of various approximations on the outcome of 
the simulations. The models can be directly compared to Lagrangian models of 
Murray (1996, 1998). We perform both isothermal simulations, and simulations 
in which the full energy equation with a realistic cooling term is solved. We 
also compare "rainfall- type" mass transfer models with those based on realistic 
modelling of the stream from the secondary. The physical assumptions on which 
our models are based are described in Sect. 2 together with numerical methods 
emplyed to solve the equations of hydrodynamics. The models are presented in 
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Sect. 3, and the results are discussed in Sect. 5. 



2 Physical assumptions and numerical methods 

We simulate the flow of gas in the orbital plane of a binary consisting of two 
stars in circular orbits around the center of mass. We use spherical coordi- 
nates end a corotating reference frame centered on the primary, with the z-axis 
perpendicular to the orbital plane. Assuming that the ratio H/r of the disk 
is constant, and the latitudinal velocity component is negligible, we can write 
the continuity equation and the equations describing conservation of radial and 
angular momentum in the following form: 

|+V.(pv) = 

where for any variable a 



V • (av) = 



1 dr'^avr 1 
+ - 

^2 Qj. J, 



In these equations j = rpV(i) is the angular momentum density measured in 
corotating frame, and f is the external force (gravitational and inertial) acting 
on unit volume. Viscous forces F^*^*^ are given by standard formulae (see eg. 
Landau and Lifshitz 1982). We assumed that kinematic viscosity coefficient and 
bulk to shear coefficients ratio are constant troughout the disk. 

The models are either isothermal or radiative. In the first case we use an 
isothermal euation of state 

P = cIp, 

where Cg is the isothermal sound speed, assumed to be constant in space and 
time. In the second case the equation of state of an ideal gas with the ratio of 
specific heats 7 equal to 5/3 is used, and the energy equation 

3F 

— + V . (Ev) = -pV . V + Q^''^ - Q'^^ (1) 
at 

is additionally solved, with E standing for the internal energy density, Q^*^^ - 
for heat generated by viscous forces and Q^"^ - for heat radiated away. Cooling 
processes are described in the same way as in Rozyczla and Spruit (1993). 

The equations are solved with the help of an explicit Eulerian code described 
by Rozyczka (1985) and Rozyczka and Spruit (1993). The inner edge of the grid 
is located at r^n = 0.1(i from the primary (where d is the distance between the 
components of the system), and the outer one {Rout) - at the inner Lagrangian 
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point Li. For all models the same mass ratio q = 3/17 is assumed, resulting 
in r^i = 0.736d. The spacing of the grid is logarithmic in r and uniform in 
the azimuthal coordinate 0, and there are 50 and 60 grid points in r and 0, 
respectively. 

The gas can freely flow into the computational domain or out of it through 
the outer boundary of the grid. At the inner grid boundary only free outflow is 
allowed for, and no inflow can occur. Due to the explicit character of the nu- 
merical scheme, the length of the time step is limited by the Courant condition. 
In the present simulations the Courant factor is 0.3. 

3 Initial setup and results of simulations 

In the present paper we repeat and refine two simulations performed by Murray 
(1996), also referred to as models 4 and 5 in a later discussion by Murray 
(1998). As mentioned in the Introduction, Murray uses a Lagrangian SPH 
code, whwreas we work with an Eulerian, grid based one. Additionally, we 
employ a more realistic physical scenario, taking into account the effects of 
cooling, and using an ideal gas equation of state instead of an isothermal one. 
Following Murray we adopt g = 3/17, a kinematic viscosity coefficient of 2.5 x 
10~^(PQ^ (where is the orbital frequency of the binary), and a sound velocity 
of 0.02 d for the isothermal cases. 

As an initial condition we chose isothermal, Keplerian disk with density 
varying as r~^-^. The disk extends from the inner boundary of the grid up 
r = O.Sbd. The rest of the grid is originally filled with a rarefied ambient 
medium at Keplerian rotation around the primary. The density of the ambient 
medium is 1% of the disk density at r = 0.35, and its pressure matches that of 
the disk. 

The original two models of Murray differ only in details of mass transfer. In 
the first case the gas from the secondary enters the computational domain as a 
narrow stream originating at Li; in the second case the gas "rains" uniformly 
onto the orbit corresponding to the circularization radius. In the following, 
they are referred to as stream and rain models, respectively (note that, by the 
definition of the circularization radius, if the mass transfer rate is the same in 
both cases then the rate of angular momentum transfer is also the same) . Below 
we describe four models: radiative stream, radiative rain, isothermal stream, 
and isothermal rain. Hereafter, they are referred to as Models (or Cases) 1-4, 
correspondingly. 

All models were followed for 100 orbits of the binary. In order to check if 
the balance between the matter added to the disk and accreted onto the white 
dwarf is established in the course of evolution, we monitored the disk mass as 
a function of time (Fig.l). For models 1, 2 and 3 the mass of the disk grows 
rapidly in the beginning, and then it reaches an equilibrium value. That value 
is almost two times greater in Case 4 (isothermal rain) than in Case 2 (radiative 
rain). In Case 3, after the phase of initial growth, the mass of the disk begins to 
oscillate with a period of ~ 45P (where P is the orbital period of the binary). 
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Figure 1: The mass of the gas contained in the grid as a function of time. The plots 
are normahzed to the initial mass of the disk. 



The next step was to check, whether the disks in our models became eccen- 
tric. For natural numbers k,l we can calculate 



/rout 
Pk,i r'^dr 



where pk^i are defined by equation 



Pk,i exp[i(/c6> - IVt^t)] 



and the term results from the assumption of constant ratio H/r (in practice, 
the lower integration limit can be equal to r^n, whereas the upper one has to 
be placed at r = 0.6 in order to avoid highly asymmetric contribution from 
the stream). If the mass of the disk is constant, then an appropriate measure 
of the eccentricity is Si^q. If, however, the mass is varying in the course of 

def 

simulation, S[ q should be used instead, where S[q = Si^q/ Sq^q (note that 5*0,0 
is proportional to the mass of gas in the integration domain). Fig. || shows 
S'lQ di function of time. It can be immediately seen that S[ o is much greater 
in isothermal models than in radiative ones. In other words, isothermal disks 
more willingly become elliptic. Secondly, only Model 4 (i.e. isothermal rain) 
exhibits a phase of exponential growth of S[ o which was predicted by Lubow's 
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Fi gure 2: 5'i,o/*S'o,o as a function of time. 



theory. Even in that case, however, another prediction of the theory is not met, 
namely, the behavior of the (2,3) mode. The strength of that mode should be 
proportional to the time- derivative of the strength of (1,0) mode, whereas in our 
simulations it is roughly proportional to the strength of (1,0) mode (see Fig. ||). 

Following Murray (1996, 1998), we calculated the rate of viscous dissipation 
in the disk as a function of time. The results of Foureir analysis of that function 
for Models 3 and 4 in range of 60 — lOOP are shown in Fig. ^ and ||. In Case 
4 (isothermal rain) the rate of viscous dissipation reaches a stationary value, 
whereas in Case 3 (isothermal stream) it exhibits large amplitude oscillations (by 
a factor of 2) with a period of ~ 45P. Note that the maxima of the dissipation 
rate coincide with the maxima of the the disk mass. In both models much 
more rapid oscillations of the dissipation rate are also seen, whose dominant 
frequency is a little lower than 1P~^ (all remaining peaks in power spectra can 
be interpreted as its higher harmonics). The corressponding periods are equal 
to l.llP in Model 3., and 1.08P in Model 4. In Model 3 the high-frequency 
oscillations are visible only within dissipation rate maxima (on the ascending 
branch shortly before peak, and on the entire descending branch). In Model 4 
they are excited when the mass of the disk approaches its maximum, and they 
persist till the end of simulation. A detailed examination of i high-frequency 
oscillation reveals their similarity to those found by Murray (1998) in his Models 
4 and 5 (see his Fig. 6). Because they are visible both in model with and without 
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Fi gure 3: 5'2,3/*S'o,o as a function of time. 



stream, those oscillations cannot be excited by the interaction of the stream with 
the advancing and retreating edge of the disc, and they most probably must be 
related to the tidal influence of the secondary. We found no similar oscillations 
in our radiative models. 

In Fig. H and ^ we present sequences of disk density maps for Models 1 
and 3 (radiative stream and isothermal stream), covering the time-span from 
99P to lOOP (at that time the rapid oscillations of viscous dissipation rate are 
fully developed in Model 3). In the radiative case the disk is almost circular. 
However, in the isothermal case the disk has a clearly nonaxisymmetric shape, 
which, at least in some frames, does not significantly differ from the elliptical 
one. The ellipse is precessing with a period a little longer than IP. In one 
precession period the maximum viscous dissipation occurs when the major axis 
of the ellipse is perpendicular to the line connecting the components of the 
binary (t = 99.65P in Fig. 0) (the same results were obtained by Murray). 
In both cases a nonaxisymmetric feature in the form of two spiral arms is also 
seen, which in Case 1 remains stationary in the corotating frame. We interpret 
it as the (2, 2) mode, excited and maintained by the tidal interaction of the 
secondary. 

Fig. ^ shows the density map of Model 3 at T = SOP (i.e. between the 
large-scale maxima of the dissipation rate, when no short-period oscillations are 
visible). At that moment the disk is nearly axisymmetric, and only the (2,2) 
mode has a significant amplitude. Thus, we conclude that the overall behaviour 
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Figure 4: Viscous dissipation rate (program units) in Model 3 as a function of time 
(top), and its power spectrum in time range of 60 — lOOP (bottom), dissipation rate 
in this time range. 
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Figure 5: Viscous dissipation rate (program units) in Model 4 as a function of time 
(top), and its power spectrum in time range of 60 — lOOP (bottom). 
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T=99.8P T=100.P 



Figure 6: Density maps for Model 1. The snapshots are taken at t = 99. P, 99. 2P, 
99.4P, 99.6P, 99.8P and lOOP. The white dwarf and the Li point are at the center, 
and in the middle of the lower edge of each frame. 
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T=99.8P T=100.P 



Figure 7: Density maps for Model 3. The snapshots are taken at t = 99. P, 99. 2P, 
99.4P, 99.6P, 99.8P and lOOP. The white dwarf and the Li point are at the center, 
and in the middle of the lower edge of each frame. 
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of Model 3 is in agreement with the tidal instability model of superhumps. 

4 Discussion 

In the present paper we reported a series of 4 simulations of disk in Cata- 
clysmic Variables of SU UMa type. The least realistic model was obtained with 
an isothermal equation of state, and a " rainfall- type" approximation for mass 
transfer from the secondary. The most realistic simulation involved a full energy 
equation, and a stream of matter originating at the inner Lagrangian point. 

We found that only the isothermal disks exhibit a clear tendency toward el- 
liptic distortions accompanied by precession. Within the framework of Lubow's 
(1992) theory, one could try to explain this result by assuming, that the strength 
of the tidally excited (2,2) mode is significantly larger in radiative models than 
in the isothermal ones. This is because, according to the theory, the (2,2) oscil- 
lations tend to keep the disk gas away from the 3:1 resonance responsible for the 
growth of the tidal instability. Such an assumption, however, is not confirmed 
by the analysis of our results: the (2,2) mode appears to be equally strong in 
all models. 

The phase of an exponential growth of the (1,0) mode, foreseen by the theory, 
was found only in the least realistic model. However, contrary to theoretical 
predictions, even in that case the time derivative of the strength of (1,0) mode 
was not proportional to the strength of the (2,3) mode. 

Oscillations with a period slightly longer than P, which may be tentatively 
associated with superhumps, were observed in isothermal models only. The 
period of oscillations found in radiative models was about 3 times shorter. 

The main results of this work may be summarized in three points: 

1. As foressen by the tidal instability theory, isothermal develop an appre- 
ciable eccentricity, and begin to precess. The precession period is the same as 
the period of rapid fluctuations in viscous dissipation rate, and it is slightly 
longer than the orbital period of the binary. The behaviour of the (1,0) mode, 
however, is not consistent with the theory. 

2. In radiative models an elliptic, precessing disk does not develop. The 
dominant oscillations in the viscous dissipation curve have periods of ^ 1/3P. 

3. In all models the two-armed (2,2) mode, excited and maintained by tidal 
forces of the secondary, is very clearly seen. 

Our conclusion is that the mechanism of superhump phenomenon is not yet 
entirely understood, and further research on this cubject is desirable. Rozyczka 
& Spruit (1993) simulated a viscosity-free disk with radiative losses, in which 
angular momentum was transported by spiral shocks. They found an erup- 
tive instability, qualitatively similar to outbursts of Dwarf Novae. During the 
eruption, the radius of the disk increased substantially. We suggest that low- 
viscosity disks prone to this type of instability may develop eccentricity and 
exhibit superhump- like oscillations during outbursts. This issue is presently 
under investigation. 
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